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Abstract 

The time evolution of a many-fermion system can be described 
by a Green's function corresponding to an effective potential, which 
takes anti-symmetrization of the wave function into account, called 
the Pauli-potential. We show that this idea can be combined with 
the Green's Function Monte Carlo method to accurately simulate a 
system of many non-relativistic fermions. The method is illustrated 
by the example of systems of several (2-9) fermions in a square well. 



1 Introduction 

The application of Green's Function Monte-Carlo (GFMC) algorithms for the 
simulation of bosonic and fermionic systems is well known , @ . However, 
the fermionic case is much more difficult to deal with than the bosonic one [Q] , 
[0], 0. In the framework of simulations of many-fermion systems employing 
the Langevin equation, Tursunov and Zhirov ^ introduced the idea of a 
Pauli-potential, to account for the repulsive forces between fermions due to 
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anti-symmetrization. We study the implementation of this idea in the more 
efficient GFMC method. 

In this work we ffist describe the standard GFMC method, then we dis- 
cuss the proposal of Tursunov and Zhirov. After that we illustrate the im- 
plementation of the proposed approach within the frame work of GFMC 
calculations of the simple quantum mechanical system of several fermions in 
a square well. With moderate computational effort we simulate nine spinless 
fermions 

2 The Green's Function Monte Carlo Algo- 
rithm 

In the paper p the authors applied the Langevin-equation to study multi- 
fermion systems. We use a more efficient method: Green's Function Monte 
Carlo 1^, in a modified form ||^. The idea of the GFMC- method is to 
determine the ground-state energy and wave function by operating iteratively 
with the Green's function on an arbitrary function. The Green's function 
itself is obtained as a solution of the standard resolvent integral equation: 

G{E) = Gt{E) + G{E){V ~Vt)Gt{E). (1) 

The (exactly known) Green's function Gt{E) is the resolvent for the Hamil- 
tonian with the potential Vr- The time evolution of the system is determined 
by the propagator in the time representation. For imaginary time this propa- 
gator is the Laplace-transform of G{E). Because we solve eq. (|I|) by iteration 
using the standard MC method, in order to achieve fast convergence, it is 
important to employ a trial potential Vr that is as close as possible to V. Be- 
sides this trial potential, there are other elements in the GFMC-method that 
make that method so efficient compared to other stochastic methods. These 
features are: the guidance function that guides the MC process and the 
trial energy Et] the details of the standard GFMC-method are described in 
the Appendix. 

We use the modified algorithm proposed in ref.0; which allows to work 
with the integral equation of the type (|l]) even in the case where the ker- 
nel is not positive definite. In the standard approach the wave-function is 
represented by a set of points, which move randomly and may disappear 
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or reappear with some multiplicity, proportional to the kernel of the integral 
equation (see step 8 of the algorithm described in the Appendix). In the mod- 
ified GFMC method the multiplicity mj (see ( |A.8D ) is proportional to the 
absolute value of the kernel, and all points which are going to the intermediate 
generation with that multiplicity, also change their phase: S{x) S{x) + Sk, 
where 6k is the phase of the kernel K which enters into the definition of m/ 
( |A.8|) . To explain this modification we give a simple example. To calculate 
the "expectation value" < >= / (j){x)f{x)dx/ J f{x)dx for the complex 
function f{x) we can generate the set of points {xk} with the probability 
proportional to |/(a;)|, and calculate < >= J2k4>{^k)e-'^^'' / J2k^^^'' where Sk 
is the phase of f{xk). Of course, the convergence of this procedure is not 
guaranteed for all choices of and /. Our results show that for the case stud- 
ied here (fermions in a square well) there exists a sufficiently broad range of 
parameters of the algorithm, in which convergence is obtained. 

The simple improvement described above allows for inclusion of the sign 
of the wave function: the density of points corresponds to the magnitude 
of the wave function. The phase of a point corresponds to the phase of 
the wave function at that particular position. In this way we are also able to 
accommodate dynamical nodes in the wave function, which occur in fermionic 
systems already in the ground state. Our algorithm provides the dynamical 
nodes in the wave function owing to the action of the Pauli-potential (|^). 

3 The Idea of Tursunov and Zhirov 

The main complication for the use of stochastic simulations of multi-fermion 
systems is the fact that its wave function must be ant i- symmetrized. A basic 
tool of the methods used here is the imaginary-time single-particle propagator 

/ r ■ \ f m fx-^ — x™) /x-^-l-x*"\l 
U (x/, X-; P)=C exp " j ' 

where m is the mass of the particle and /3 is a (small) time step. For two 
identical fermions, 1,2, the imaginary-time propagator can be written as: 

[/(^)(x(,x^,xr,x-;/?) = 
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(3) 



where (F) refers to fermions and (D) to distinguishable particles. U^^'^ is the 
product of two single-particle propagators: 



f/(^) (xf , x^, xf, x^"; P)=U (x{, xi"; /3) U {4, x: 



in 
2 1 



(4) 



At small values of f3, anti-symmetrization (^) can be effectively implemented 
[|] by using an additional effective potential: 

f/^(l,2) = f/^(l,2)-f/^(2,l) = f/^(l,2)(l-f/^(2,l)/f/^(l,2)) 

^ ?7^(l,2)e-'3^"(i'2). (5) 

This additional potential, V^^\ has the following form to leading order in /5: 



V^(^)(l,2) = --ln 



1 — exp 



m(x{ — X2) 



(6) 



We see that the Pauli-exclusion principle leads to a complex, nonlocal and 
time-dependent potential; still, it can be used in the Monte-Carlo algorithm. 
For the N-fermion case we will have 



V^^\1,...,N) 



(7) 



— In 



1 - exp - Efc 



So, to leading order in P, this "Pauli-potential" corresponds to anti- 
symmetrization of pairs of particles only: the sum in eq. (|^ has only N(N- 
l)/2 terms. It is very important that the permutations of three and more 
particles occur when the time-development of the system is simulated by re- 
peated operation of the propagator U^^\ i.e., by the repeated action of the 
potential p during the Monte Carlo procedure. 



4 Results 

We obtained our results by solving eq. (|^) in the time-representation. Using 
this Green's function, we obtain the ground-state wave function of the system 
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dependent on time: ipoix; [3) ~ exp(— i?o/?)V'o(a;; 0). For an N-particle system 
in three dimensions, the wave function is represented in the MC-method by 
a set of points in 3N-dimensional space. The density of the points in such 
a population follows modulus of the wave function \'iIjq{x] [3)\. The action of 
the density matrix (Laplace-transformed Green's function) is performed in 
finite imaginary time steps (3. 

In order to obtain the Green's function, we perform a Laplace transform 
of the density matrix by sampling (3 from the distribution exp(— /3 /A) /A. We 
varied A in order to be able to extrapolate our results to the point A = 0, 
which means that, on average, the time step j3 tends to zero, since we used 
the Pauli-potential only to leading order in (3. 

In order to determine the energy Eq of the ground state, we monitor the 
size of the population in time and use ( [A.13|) from the Appendix. 

We found that we could decrease the fluctuations in the energy by an 
order of magnitude if we would kill all points that have a multiplicity (see 
items 4 and 5 of the GFMC algorithm described in the Appendix) larger 
than some value Mmax (Typically Mmax ~ 5 — 20). By this procedure we 
introduce a systematic error proportional to A. Doing so, we do not change 
the character of the systematic error: it remains linear in A. 

For a check of our approach we have, until now, studied the case of several 
spinless fermions in a square-well potential. The number of particles in this 
well was varied between two and nine. 

For the trial potential Vt we use the oscillator interaction, for which the 
Green's function is known in closed form (|A.6|) [0 . For the guidance function 
ipG we use a Slater determinant of harmonic-oscillator wave functions. (These 
may or may not correspond to the same oscillator as is used for the trial 
potential.) 

In Fig.l we show for a value of A = 0.0005 the development of the 
energy for a system of nine fermions with the number of time steps. In this 
calculation we use a square well with depth Vq = —3.5, radius R = 2; the 
number of points in each population (which describes the nine fermion wave 
function) was approximately one thousand. Clearly, the energy converges to 
a value of -12.8, which differs from the true value (-11.501), indicated by the 
broken line. The reason for this phenomenon is the fact that we still have a 
finite A. 

In order to see the effect of taking smaller time steps, we calculated the 
energy for different values of A. Fig. 2 shows this dependence of the average 
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energy for the same problem on the size of A. We clearly see that the average 
energy tends to the exact energy if A tends to zero. Extrapolation of the 
energy values to the point A = gives a value that is almost equal, within 
the error bars to the exact value for this case (-11.501). These simulations 
were performed for a value of Ex = —10.5, differing from the exact value of 
the energy, to mimic a realistic situation, when the value of the exact energy 
is not known. 

The true value of the energy can be found by extrapolating the computed 
values to A = 0. We illustrate this for the case of nine bodies in Fig. 3. The 
number of killed points (taking into account their multiplicities) divided by 
the total number of points in our simulation process is plotted as a function 
of A. We see that this number indeed depends linearly on A. A value 
Mmax = 5 was taken. We checked that a similar linear dependence occurs 
when we change M^ax- 

The dependence of the average energy on A is shown in Fig. 4 for the case 
of five fermions in the same square-well potential. The linear dependence of 
the energy on A is clearly seen to occur for sufficiently small values of A. 

We also found that the Pauli potential acts in such a way that the motion 
of points that render the argument of the logarithm in eqs. (^0) negative, is 
hindered. We checked, by recording the number of points that do cross the 
border, that the Pauli repulsion effectively blocks the crossing, for A(3 — >■ 0. 
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Appendix A 



Below we describe the main steps of the algorithm 0] (see also the pedagog- 
ical paper |^ or the book by Kalos and Whitlock |^ ) for the solution of the 
integral equation (|l|). 

1. A set of points is sampled from a distribution = |\E'G(a;)P- For 
fermions in D dimensions X IS cL DA^-dimensional vector representing 
the positions of A^ particles (if we exclude the center of mass motion 
it is an (A^ — l)D-dimensional vector). For the initial distribution we 
take |\E'G(a:)P, because eventually we will obtain from the process the 
distribution ^'(x)^g(x) rather than \l/(x). The typical number of points 
in this set, called the initial generation, is several hundreds. 

2. It is convenient to work with the density matrix p(x, = 
J2n'^n{^)'^n{x')e~^"^ . The density matrix is related to the Green's 
function G{x, x', E) by the Laplace transform. In our algorithm we 
carried out the Laplace transform by sampling the imaginary time (3 
from the distribution exp | ~ a}' ^u^d that the additional random 
number occurring in the sampling of the imaginary-time distribution 
involved in the Laplace transform, improves the statistical errors. 

3. To each point x' in the initial generation the diffusion and drift is 
applied, after which the points are distributed with the probability 

f{x,f3) = J pn{x,x',P)Mx')dx', (A.l) 

where 

DN 

, ^. ( m \ { mix — x' — BF{x')Y\ , 

m is the mass of the fermion and the "quantum force" F is given by: 

F,(x) = ^'^g\^)^^g{x), (A.3) 
Note that Fi{x) is the component of a vector cls is. 
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Because po does not enter the final expression, as the dependence on 
Pd is cancelled due to the fact that pd enters also in the denominator 
of rriD (see ( |A.4| )) the random walk of points can be performed with an 
arbitrary probability distribution. For the distance probability distri- 
bution a Gaussian form is used, since random numbers with a Gaussian 
distribution can be generated very efficiently on a computer. The sec- 
ond reason for choosing this distribution is, that it is similar to the real 
distribution corresponding to the bound state. 

Algorithmically the sampling of the distribution (|A.1|) is done in two 
steps. First we shift the initial points: x" = x' + (3F{x'), and after 
that we add to each point the gaussian random numbers r]i with unit 
expectation value: Xi = x" i + \f^rii 

4. In order to construct the new generation, multiple copies of each point 
are produced. The multiplicity is given by the formula: 



X 



where pt is the density matrix (the Laplace-transformed Green's func- 
tion) for the trial potential which satisfies the well known equation: 

^PI^^^^ = -Htp{x,x',P)- pt{x,x',0)=6{x-x'). (A.5) 

For the trial hamiltonian we use the harmonic oscillator hamiltonian, 
for which px is known explicitly 0: 



X exp < ^'^T r/ 2 _|_ ^/2\ QQgj-^^y^ _ 2x ■ x] — ct(3 \ 

[ 2 smh UtP J 

This trial density matrix corresponds to a hamiltonian with the poten- 
tial Vt, given by 
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Vt{x) = + ct. (A.7) 

The trial energy Et is introduced in order to avoid exponential growth 
or shrinkage of the number of points of the population. As a result, the 
number of points in the population fluctuates around a value that has 
the time-dependence exp{ — (£^0 — Et)[3}. Effectively, this amounts to 
a shift in the hamiltonian with the constant energy Ex- 

5. The other, so called intermediate, branch of the process is formed by 
creating another set of multiple copies of the points x: for them the 
multiplicity is given by: 

, K{x,x\(3)mD{x,x',l3)/^ 

mi{x, X ,(3) = — , A.8 

Pt{x,x',I3) 

here K is the Laplace transformed kernel of the integral equation (|l]). 



K{x, x', 13) = [Vrix] - V{x)]pTix, x', /3). (A.9) 

In general mo and mi are not integers. We convert them to integers 
by adding a uniformly distributed random number to each of them and 
take the integer part. 

6. Each intermediate point, created this way, is treated in the same way as 
the points taken from the initial generation, i.e., they will take part in 
the random walk with branching until they are eventually propagated 
to the new generation. If the average value of mj is less than unity, 
this process is completed in a finite time on the computer. 

The sequence of operations described above, correspond to the terms in 
the iterative solution of the equation: 

P = Pt + A-K*p (A.IO) 

where K = V and * denotes the integration over the intermediate coordi- 
nates: K * p = J K{x,x")p{x" ,x') dx" . The direct points correspond to the 
term pr, while the intermediate points correspond to A ■ K. If an interme- 
diate point is processed again, it may be promoted immediately to the new 
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generation, in which case it corresponds to the term A ■ px * K, otherwise it 
will be an intermediate point again, now corresponding to ■ * etc. 
The difference in the normalization of the kernels of the equations (|I|) and 
([A. 10|) is due to the difference in the normalization of G and p: 



here \l/n(a^) and En are the wave function and the energy of the n-th level of 
the Hamiltonian studied here. The distribution of points in the new (second) 
generation is sampled with the probability distribution 

/2(x) = / p{x,x')fi{x')dx'. (A.12) 



In the same way we can get the distributions /s, /„; 

Using expression ( |A.11| ) for p{x, x') it is easy to prove that as n ^ cxo 
fn{x) const. (a;) ^0 (a;) /[I + A(Eo - ^t)]"~^ + where ^^0(2;) and Eq 
are the wave function and the energy of the ground state. So the ground 
state energy can be calculated from the number of points, P„ in the n-th 
generation: 

E, = Et + ^(^-i]. (A.13) 



A V Pn 
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Figure Captions 

Fig. 1 The development with the number of time steps (N) of the energy for 
a system of nine fermions. The value of A is 0.0005. 

Fig. 2 Dependence of the energy of the system of nine fermions on the time 
step A. A linear fit is made, which gives an estimate of the energy for 
A = 0. 

Fig. 3 The dependence on A of the number of killed points for the system of 
nine fermions divided by the total number of points in our simulation. 
The parameter M^ax — 5. 

Fig. 4 The same as Fig. 2, but now for five fermions. It is seen that the 
linear dependence of the energy on A obtains for small values of A. 
A linear fit of the points at low values of A intersects the axis at the 
exact energy (indicated by the broken line). 
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